######################################################
# CBQ functionality
######################################################

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

pald <- function(x,p){
  out = ifelse(x<0,p*exp(x*(1-p)), 1-(1-p)*exp(-x*(p)))
  return(out)
}

cbq <- function(N, # number of obs
               n_covariate, # number of covariates
               X, # covariates
               Y, # d.v
               N_indx, # number of choice sets
               ind, # index of choice sets
               qtl, # quantile to be estimated 1...9
               nchain = 2, # number of parallel chains
               niter = 1000, # number of interations
               thin = 1, # thinning
               seed = 1234567, # random seeds
               offset = 1e-20 # enhance numerical stability
               ){
  datlist = list(N=N,
                 Y = Y,
                 D_common = n_covariate,
                 X_common = X,
                 N_indx = N_indx,
                 ind = ind,
                 offset = offset)
  pars = c("beta_common")
  init = function(){
    list(beta_common=rnorm(n_covariate,0,1)
    )
  }
  ms_cbq_model = stan_model(paste("replication_MS/stan_new_models/cbq4.0_q",qtl,".stan",sep=""))
  ms_cbq_stan = sampling(ms_cbq_model, 
                         data=datlist,
                         pars=pars, 
                         seed=seed,
                         chains=nchain,
                         iter=niter,
                         thin=thin,
                         init=init)
  return(ms_cbq_stan)
}
